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Abstract 

In this paper we consider the classic problems of scattering of waves 
from perfectly conducting cylinders with piecewise smooth bound- 
aries. The scattering problems are formulated as integral equations 
and solved using a Nystrom scheme where the corners of the cylinders 
are efficiently handled by a method referred to as Recursively Com- 
pressed Inverse Preconditioning (RCIP). This method has been very 
successful in treating static problems in non-smooth domains and the 
present paper shows that it works equally well for the Helmholtz equa- 
tion. In the numerical examples we specialize to scattering of E- and 
H-waves from a cylinder with one corner. Even at a size kd = 1000, 
where k is the wavenumber and d the diameter, the scheme produces 
at least 13 digits of accuracy in the electric and magnetic fields every- 
where outside the cylinder. 

1 Introduction 

The numerical simulation of scattering from cylinders has a long history in 
computational electromagnetics. As early as in 1881, Lord Rayleigh treated 
the scattering of light from a circular dielectric cylinder [1] . He considered an 
incident plane E-wave, i.e., the electric field is parallel to the cylinder, and 
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a permittivity and permeability of the cylinder that departed only slightly 
from those of the surrounding medium. That enabled him to find an approx- 
imate solution that today is referred to as the Born approximation and can 
be viewed as spectral method solution with only one basis function, c.f. [21 
Section 8.3.4]. The theory of scattering from circular cylinders and spheres, 
conducting or dielectric, was soon after that fully understood by using expan- 
sions of the incident and scattered waves in partial waves, c.f. [3]. Since then, 
a large number of papers have been published that solve scattering problems 
in electromagnetics, as well as in acoustics and elastodynamics, using differ- 
ent numerical techniques. All with the common goal of constructing faster 
and more accurate solvers for ever more detailed and complex geometries in 
two and three space dimensions. In particular, integral equation methods 
have become very important tools. In electromagnetics such methods were 
made popular by the contributions of Harrington, c.f. [4]. The mathematical 
foundations of the scattering problems and the integral equation formulations 
are discussed in the books by Colton and Kress [6] . 

The present paper is about scattering from piecewise smooth perfectly 
conducting objects. The presence of boundary singularities, such as corners, 
tends to cause complicated asymptotics in quantities used to represent the 
solution. Intense mesh refinement might be needed for resolution, but this is 
costly and can easily lead to instabilities and the loss of precision in the com- 
puted field. In the context of integral equation solvers, regions close to the 
boundary are the most problematic. On the application side, scattering from 
non-smooth metal objects is of great importance in radar imaging of objects 
with sharp corners such as airplanes, vessels and vehicles. Sharp corners 
that are oriented perpendicular to the line of sight of a monostatic radar 
may create reflections that are large enough to be detected by the radar. 
The two-dimensional approximations can be used for elongated objects like 
wings but also in the evaluation of fields in the near zone of smaller objects. 
Other important two-dimensional problems are wave propagation in rectan- 
gular waveguides, photonic band gap structures, and substrate integrated 
waveguides. 

The numerical solver used in this paper takes its starting point in a Fred- 
holm second kind integral equation with integral operators that are compact 
away from boundary singularities and whose unknown quantity is a layer den- 
sity representing the solution to the original problem. The integral equation is 
discretized using a Nystrom scheme and composite Gauss-Legendre quadra- 
ture. At the heart of the solver lies a method called Recursively Compressed 
Inverse Preconditioning (RCIP). It modifies the kernels of the integral opera- 
tors so that the layer density becomes piecewise smooth and simple to resolve 
by polynomials. Loosely speaking one can say that RCIP makes it possible 
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to solve elliptic boundary value problems in piecewise smooth domains as 
cheaply and accurately as they can be solved in smooth domains. The RCIP 
method originated in 2008 [7J. In a series of papers it has been extended and 
successfully applied to electrostatic and elastostatic problems which, at first 
glance, might seem outright impossible. For example, the effective conduc- 
tivity of a high-contrast conducting checkerboard with a million randomly 
placed squares in the unit cell was computed on a regular workstation with 
a relative accuracy of 10 -9 [8]. A new record has been established for the 
three-dimensional problem of determining the capacitance of the unit cube 
- 13 digits compared to the seven digits that were previously known [9]. 

When we here apply the RCIP method to the Helmholtz equation and the 
exterior Dirichlet and Neumann problems we do this in a two-dimensional 
setting. We consider scattering of time-harmonic E- and H-waves from an 
infinitely long perfectly conducting cylinder. Scattering problems are harder 
to solve than electrostatic problems, all other things held equal. Planar 
problems provide a good testing ground prior to a move up to three dimen- 
sions [10]. As we shall see, the transition from Laplace's equation to the 
Helmholtz equation is surprisingly straightforward and the results, presented 
in Section [4] below, are as good as the ones obtained for electrostatics. 

Our numerical solver meets five important criteria. The first criteria is 
that it can handle cylinders with general shapes. In practice this means 
cylinders with piecewise smooth boundaries and with a finite, but arbitrary, 
number of corners. The second criteria is that it can treat frequencies ranging 
from zero up to large values of kd, where k is the wavenumber and d the diam- 
eter of the object. We have found that kd = 1000 is quite easy to reach and 
for most cylinders this frequency range overlaps the frequency band where 
approximate high frequency methods, e.g. unified theory of diffraction in 
combination with physical optics, can be applied with reasonable accuracy. 
The third criteria is that the method can deliver accurate results for the 
scattered field everywhere outside the object. Even close to a corner and at 
kd = 1000 the scattered field is calculated with at least 13 digits of accuracy 
in IEEE double precision arithmetic (16 digit precision). The fourth criteria 
is that the method enables fast solvers. In the present implementation the 
solver is fast only in the sense that the cost for modifying the kernels of the 
integral operators grows linearly with the number of corners in the compu- 
tational domain. The method can be made fast in toto by incorporating fast 
multipole techniques [UJ, [12] or perhaps even fast direct solvers [13], E] . The 
fifth criteria is that the method is automatized and flexible. It requires only 
a minimum of adjustments as operators and geometries change. 

It is beyond the scope of the present paper to review the RCIP method 
in its entirety. In Section [3] we give a brief overview and a few details on 
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discretization issues particular to Hankel kernels. Apart from that, we refer 
readers to the original research papers [TJ [151 EH IZ] and to a newly written 
tutorial [18J. 

There are several recent journal papers that focus on speed and accuracy 
for two-dimensional scattering problems in complex geometries. In [19] scat- 
tering from two-dimensional smooth strips are treated using integral equa- 
tions and a Nystrom method. In [20] the approach of [19] is generalized to 
smooth slotted cylinders. A similar problem is treated in |21j . The schemes 
used in these papers give accurate results but they cannot, in a simple way, be 
generalized to non-smooth geometries. In [14] and in [22J, on the other hand, 
very fast and also flexible and accurate numerical schemes are developed for 
the solution of integral equations modeling scattering from general objects 
with both corners and multi-material junctions. These papers, however, do 
not address the problem of accurate near field evaluation. 

2 Formulation of the problems 

We consider in-plane waves scattered by a bounded perfectly conducting 
cylinder with a piecewise smooth boundary T. The region outside the object 
is denoted fl eyL , the time dependence is e~ luJt and r = (x,y). Both E- waves, 
often referred to as TM-waves, and H-waves, often referred to as TE-waves, 
are treated. We decompose the electric and magnetic fields into a sum of 
the incident field, denoted U inc (r), generated by a source in f2 ex , and the 
scattered field, denoted U sca (r) in both cases. 

2.1 E- waves 

We let the electric field be parallel to the cylinder, E(r) = zll(r), and let 
U{r) = U mc {r) + U SCSL (r). The scattered field t/ sca (r) satisfies the following 
exterior Dirichlet problem: 



We write the solution as the combined integral representation [6l eq. (3.25)]. 



V 2 f/ sca (r) + k 2 U sc& {r) = 0, r E 
U s Ur) = -U- mc (r), r ET 



■ex 



(1) 

(2) 




(3) 




r Gfi, 



(4) 
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where $ k (r,r') = -iJg {k\r — r'\) is the free space Green function for the 

Helmholz equation in two dimensions, Hq is the Hankel function of the first 
kind of order zero, and d£ is an element of arc length. The index k indicates 
that the quantity or function depends on the wavenumber k = u/c. Insertion 
of © into @ gives the integral equation for the layer density p(r) 

(I + K k - i|s fc )p(r) = -2U inc (r), r e T, (5) 

where 

K k p(r) = 2 / d ®f y) p{r')<M! (6) 
S kP {r) = 2 J <$> k {r,r')p{r')dl'. (7) 

k 

The second term on the right hand side in Q corresponds to the term l ~^S k 

in ([5]) and is added in order to ensure a unique solution for all k. The 
equation (|5]) is often referred to as an indirect combined field integral equation 
(ICFIE). 



2.2 H- waves 

We let the magnetic field be parallel to the cylinder, H{r) = zll(r), and let 
U(r) = U mc (r) + U SCSi (r). The scattered field U sceL (r) satisfies the following 
exterior Neumann problem 



V^ sca (r) + k z U sc& {r) =0,reQ 
dUarJr) dU- mr (r 



reV (9) 
lim [A- ik] U sca (r) = 0, (10) 



dv r dv, 
.im 

\r\-^oo \Or 

where — ^ ca ^ - is the normal derivative of U sca . There are several ways to 



dv. 



r 



model this problem as an integral equation. We use a regularized combined 
field integral equation since it is always uniquely solvable. The scattered field 
is then obtained from the representation 



U sca (r) = J Hr,r')p(r')de' + i J S ik p(r'W, r e fi cx , (11) 



5 



which after insertion into ^ gives the integral equation 

.dU- mc (r) 



(I-K' k -iT k S ik )p(r) = 2- 



dv r 



(12) 



Here K' k is the adjoint to the double layer integral operator K k in (|6 



and 



T k p(r) 



d 



K k p(r). 



(13) 



(14) 



The equation ( 12 ) is sometimes referred to as ICFIE-R 



It is useful to observe that the hypersingular operator T k in (14) can 



be expressed as a sum of a simple operator and an operator that requires 
differentiation with respect to arc length only [24] 



T k p(r) = 2k 2 J $ fe (r, r'){u r ■ v r ,)p(r')dt' + 2^ J $*(r, r 



df 



We may then rewrite (12) in the form 



8U- (r) 

(I + A k - iB k S ik - iC k C ik )p(r) = 2 ^ ncl \ r e T, 



dv r 



(15) 



where A k = —K', and 



B k p(r) = 2k 2 J^ k (r,r')(u r ■ u r ,)p{r)de (16) 
C k p(r)=2^J^ k (r,r')p(r>)d£>. (17) 



3 Numerical scheme 

This section briefly reviews the RCIP method, for obtaining accurate solu- 
tions to integral equations on piecewise smooth surfaces, with focus on basic 
concepts and on some details particular to the Helmholtz equation. A richer 
description, along with demo codes in Matlab, can be found in [18]. 
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3.1 Basics of the RCIP method 



Assume that we have an integral representation of a field U(r), r 6 Q ex , in 
terms of a layer density p(r) on a piecewise smooth boundary T, which leads 
to a Fredholm second kind integral equation 

(I + K)p(r)=g(r), reT. (18) 

Here / is the identity, g is a piecewise smooth right hand side, and K is 
some integral operator with kernel K(r,r') on T that is compact away from 
a finite number of corners. Let us split the kernel 

K(r,r') = K*(r,r') + K°(r,r') (19) 

in such a way that K*(r, r') is zero except for when r and r' both lie close 
to the same corner vertex. In this latter case K°(r,r') is zero. The kernel 
split (19) corresponds to an operator split 

K = K* + K°, (20) 

where K° is a compact operator. The variable substitution 

p(r) = (I + KT 1 P(r) (21) 



allows us to rewrite (18) as a right preconditioned integral equation 

p(r) + K°{I + K*)- 1 p{r) = g(r), reT, (22) 

where the composition K°(I + K*)~ l is compact. 

Let us discretize (22) using a Nystrom scheme with composite 16-point 
Gauss-Legendre quadrature. The quantities p, K°, and g should be simple 
to discretize and resolve accurately on a coarse mesh made of quadrature 
panels T p of approximately equal length. Only the inverse (/ + K*)~ x needs 
fine local meshes for its accurate resolution. We arrive at 

(I coa + K° oa R)p coa = g coa , (23) 

where the block-diagonal compressed weighted inverse matrix R is given by 

R = P^(I fin + K* n )- 1 P. (24) 



In (23) and (24) subscript "coa" indicates a grid on the coarse mesh, subscript 
"fin" indicates grids on fine local meshes, the prolongation matrix P performs 
polynomial interpolation from the coarse grid to fine grids and P^ is the 
transpose of a weighted prolongation matrix. See [18j Section 4 and 5] for 
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details. Once (23) is solved for p coa , a discrete weight-corrected version of 



the original layer density can be obtained from 

Pcoa = R-Pcoa- ( 2 5) 

The solution U (r) can then be recovered in most of the computational domain 
using p coa in a discretized version of the integral representation for U(r). 



Note that in (23), the need for resolution in corners is not visible. The 



transformed layer density p coa on a non-smooth T should be as easy to solve 



for as the original layer density p coa in a discretization of (18) on a smooth T. 
All computational difficulties are concentrated in the matrix R. Let there be 
n discretization points on the local fine grid close to a particular corner on F. 



Judging from the definition (24), it seems as if computing R should be a pro- 
hibitively expensive and also unstable undertaking for large n. Fortunately, 
R can be computed via a fast and stable recursion which relies on a hierarchy 
of small nested meshes. This fast recursion enables the computation of the 
diagonal block of R, that corresponds to a particular corner, at a cost only 
proportional to n. Actually, when very large n are needed for resolution the 
cost can be further cut down with the use of Newton's method. See jT8], 
Section 6 and 12] for details. 

The fast recursion for R can also be run backwards for the purpose of 
reconstructing p fin from p coa . A partial reconstruction of p fin is needed when 
U (r) is to be evaluated at points in f2 ex that lie close to corner vertices. 
See [TBI Section 9] for details. 

We remark that the integral equations ^ and (15), which are to be 
solved in this paper, have a more complicated appearance than the model 



equation ( 18 ). In practice this poses no problems for RCIP - just some extra 
work. The two integral operators in (|5]) can, for programming purposes, 
be combined into a single operator. The composition of integral operators 



in (15) can be treated with an expansion technique. With the help of two 



new temporary layer densities, one can arrive at a recursion for an expanded 



compressed inverse matrix R with the same structure as (24). Once R is 



computed one can extract separate blocks from it and use them in a more 



involved version of (23 ) that still uses only a single transformed global density 



p coa . See [TBI Section 14 and 17] for details. 



3.2 The discretization of Hankel kernels 

High-order accurate Nystrom discretization of boundary integral equations 
associated with the Helmholtz equation is a topic that has received much 
attention recently. See [25] for a comparison of various 2D schemes. We now 
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present our preferred scheme by showing how to discretize the operator K k 
of ^ and the first operator on the right hand side of (|4]). The other integral 
operators of Section [2] are discretized in similar ways. 

The kernel of K k is twice that of the first operator in Q and can, modulo 
a constant of i/2, be expressed as 

K k {r, r') = k\r - r'^^r - r'\) { - , (26) 

where is the Hankel function of the first kind of order one. When r G T, 
it is instructive to write (26) in the form 

K k (r, r') = f(r, r') + - log \r - r'|B {K k (r, r')} . (27) 

7T 



For a fixed r G T, we see from (26) and a series representation of that 
f(r,r') and 9ft {K k (r, r')} are smooth functions of r' G T and that 

limlog|r-r , |K{ir fc (r,r / )} = 0. (28) 
Consider now the integral / p (r) over a quadrature panel T p 

J p (r)= / K k (ry)p(r')d£'. (29) 

Let r(t) be a parameterization of T. Discretizing K k means being able to 
evaluate (29) for all r of interest, given a set of values p{r{tj)) on each T p . 

If r is a point away from T p , then K k (r, r') is a smooth function of r' G Y p 
and /p(r) can be evaluated to high accuracy using 16-point Gauss-Legendre 
quadrature 

J P ( r ) ~ Yl Kk ( r i r /)/W> ( 3 °) 

j 

where Vj = r(tj), pj = p(r(tj)), Sj = \dr(tj)/dt\, and tj and Wj are nodes 
and weights on T p . 

If r, is a discretization point close to T p or on T p , then K k (ri,r') is not 
a (sufficiently) smooth function of r' G T p and we use (27) to arrive at 

I P (n) ~ r i)Pr s i w 3 + i K k(ri, »'/!}/', (31) 

where are high-order product integration weights for the logarithmic 
operator which can be constructed using the analytic method in [TBI Section 



2.3]. The formula (31 ) can be rearranged into a particularly convenient form 

J P 0i) ~ K k( r i> r j )t>j s j"'.i + — Y ® i Kk ( r ^ T '/)}/'^/ "'/"■;}['' • ( 32 ) 

j j 
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where the weight corrections 



7 ,,corr 
W ijL 



Wjjh 
SjWj 



\og\vi - rA 



(33) 



are cheap to compute and depend only on the relative length (in parameter) 
of neighboring quadrature panels and on nodes and weights on a canonical 
panel. The formula (30) with r = r, and (32) summarize our Nystrom 
discretization of Kf. on V. 

If r is a point not on T but in fi ex close to T p , we write (26) in the form 



K k (r, r') = g(r, r') + - log |r - r'\U {K k (r, r')}+' 21 V ) ' ' V 



7T 



7T 



(34) 



We see from (26 ) and a series representation of i?{ that g(r, r') and 9ft {Kj,(r, r 
are smooth functions of r' . In analogy with (32) one can write 



2i 



7T 



Fj - r 



wj , (35) 



where u^£ rr (r) are weight corrections as in (33), but with r*j replaced by 
r, and Wjc{r) are high-order product integration weights for the Cauchy 
singular operator which can be constructed using the analytic method in [IB], 



Section 2.1]. The formulas (30) and (35) are used to discretize the first 
operator in Q when producing field plots. 



3.3 Convergence and error estimates 

Our solver shows a stable behavior. This means that the solution converges 
rapidly with coarse mesh refinement up until a point beyond which no further 
improvement occurs. Actually, beyond this optimal point there will be a slow 
decay in the quality of the solution, due to accumulated roundoff error. The 
precise location of the optimal point is hard to determine a priori. It depends 
on the geometry, on the boundary conditions, and on the wave number. The 
optimal point is determined experimentally in the numerical examples of 
Section HI 

We have estimated the accuracy in our solutions U(r) rather thoroughly. 
The tutorial [T%t Section 18] contains error plots for exterior problems in 
non-smooth domains produced in a direct way. These are achieved by gener- 
ating the boundary conditions on T via line sources inside T so that the exact 
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solution is known. In the plane- wave scattering examples of Section 4T, be- 
low, no exact results are known. Therefore we proceed as follows: we first 
compute a solution U (r) using a number of coarse panels on V deemed suf- 
ficient for resolution. Then we increase this number with 50 % and solve 
again. The difference between the resolved value of U(r) and the overre- 
solved value of U(r) is used as an indirect pointwise error estimate. Yet an 
indirect method to estimate the (overall) precision in the computations is by 
comparing the scattering cross section computed from its definition (close to 
T) with its value obtained via the optical theorem (at infinity). See, further, 
Section 4.2 As it turns out, the various error estimates seem to agree well. 



4 Numerical examples 

We shall now solve ^ and ( JT5| for the unknown density p(r), using the 
method of Section [3j and then evaluate the scattered fields of ^ and (11). 
We restrict the numerical examples to scattering from an infinite straight 
cylinder with boundary T described by 

r(t) = sin(Trt) (cos((t - 0.5)tt/2), sin((t - 0.5)tt/2)) , t G [0, 1] , (36) 

and to the incident plane wave U- mc (r) = e lky for both E- waves and H- waves. 



The object parameterized in (36) has a corner with opening angle 6 = n/2 
at r = and a diameter d — 1, in arbitrary length units, so that kd = k. 
The examples cover sizes from kd = 1 up to kd = 1000. We have seen that 
at kd = 1000 the frequency is high enough such that the uniform theory of 
diffraction theory can be applied. All numerical examples are executed in 
MATLAB on a workstation equipped with an IntelXeon E5430 CPU at 2.66 
GHz and 32 GB of memory. 



4.1 Near field 

A criterion for a powerful method is that it should be able to calculate the 
electric and magnetic fields everywhere in Q ex . Figures [T] and [2] show the total 
electric field for the E-wave and total magnetic field for the H-wave in the 
vicinity of the scattering object and the corresponding errors. The scattering 
object itself appears in green color in the left images and in white color in 
the right images. The number of spatial points in each image is 10 6 . It is 
encouraging to see, in the right images of Figures [T] and [2j that the accuracy 
is high even close to the boundary and, in particular, close to the corner. 
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Figure 1: Left: a), c), e) show dt{U(r)} for a plane E-wave L^i nc (r) = e 



iky 



incident on the perfectly conducting cylinder with boundary T given by (36) 
Right: b), d), f) show absolute errors. 
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Figure 2: Left: a), c), e) show $l{U(r)} for a plane H-wave U- mc (r) = e lky 



incident on the perfectly conducting cylinder with boundary T given by (36) 
Right: b), d), f) show absolute errors. 
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The integrals in Q and (11) are often thought of as difficult to evaluate 
close to the boundary due to the singularities in the Hankel functions when 
r' = r. However, the present method circumvents these problems using the 



high-order analytic quadrature outlined in Section |3.2| 

In Figures [l] a) , c), e) the real part of the total electric field U(r) for the 
E-wave case is plotted for kd = 10, 100, and 1000. To capture the diffraction 
pattern in the vicinity of the corner, the field is plotted in a rectangular region 
with side length proportional to 1/k and center at the tip of the corner. At 
kd = 10 the error is very small, as seen from Figure [l]b). The errors increase 
slightly with kd but even at kd = 1000 we get 14 digits or better almost 
everywhere, as depicted in Figure [I] f). For H-waves the accuracy is almost 
as good as for the E-waves, as seen from Figure [2j 

For kd = 100 and 1000 we can interpret the field plots in Figures [I] c) , e) 
and [2] c) , e) through the theory of diffraction. Thus, the outer region f2 ex is 
divided into three subregions separated by the reflection boundary and the 
shadow boundary. 

4.2 Scattering cross section and optical theorem 

In two dimensions the scattering cross section reads 

a sca = = »(- / tWO^f^df ) , (37) 

Sine -2/ Jr ciTC dv r , J 

where P sca is the scattered power per unit length, S mc -y is the component 
of the Poynting vector of the incident field, i.e. the incident power density, 
the boundary r circ is a closed curve that circumscribes the boundary T, and 
the star denotes complex conjugation. The expression holds for both E- and 



H-waves. In a numerical experiment with the cylinder of (36) we let r c i rc be 
a circle of radius 0.55 and with center at r — (0.5,0). Since the diameter 
of the scatterer is d — 1, the smallest distance between the F and r circ is 
0.05 and it occurs at the corner vertex and at a point opposite to the corner 
vertex. For evaluation points r' so close to the boundary, the field U SCSu (r') 
and its normal derivative are in general hard to evaluate. But, as we have 



already seen in Section 4.1, the RCIP method and the high-order analytic 



quadrature outlined in Section |3.2| should allow for high accuracy. 

By utilizing the optical theorem we get an alternative expression for the 
scattering cross section 



a sca = - lim^t-U^y)^^-^-^ } (38) 
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cross section E-waves 



cross section E-waves: errors 




b) 



cross section H-waves 



cross section H-waves: errors 




A) 



Figure 3: The scattering cross sections cr sca for the E-wave, a), and H-wave, 



c), calculated by the optical theorem (38) and the relative error, b) and d), 
compared to the values from equation (37) 
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which should be even simpler to evaluate than (37) since it only involves the 
far field. The mismatch between the scattering cross sections computed via 
(37) and via (38) can be used as an error estimate for both expressions. The 
cross sections for the E-waves along with such error estimates are given in 
Figure |3| a) and b) and the corresponding data for the H-waves are given in 
Figures!^ c) and d). The mismatch error is on the order of 10~ 15 . The cross 
sections in Figures [3] a) and c) show the well known behaviors for large and 
small values of k. 



5 Conclusions 

We have shown how the basic problems of scattering of E- and H-waves 
from perfectly conducting cylinders with corners can be solved numerically 
to high accuracy on a mesh that on a global level is not refined close to corner 
vertices. We give examples where the scattered electric and magnetic fields 
from a cylinder with one corner and with a diameter of up to 160 wavelengths 
is obtained with 14 digits of accuracy almost everywhere outside the cylinder. 
This success is achieved by 

1. choosing a suitable integral representation of the scattered field in terms 
of an unknown layer density 

2. formulating the scattering problem as a Fredholm second kind integral 
equation with operators that are compact away from the corners 

3. discretizing using a Nystrom scheme and a mix of composite Gauss- 
Legendre quadrature and high-order analytic product rules 

4. modifying the discretized integral equation so that the new unknown, 
a transformed layer density, is piecewise smooth 

5. solving the resulting well-conditioned linear system iteratively for the 
transformed layer density 

6. partially reconstructing the original layer density from the transformed 
layer density 

7. evaluating the scattered field from a discretization of its integral rep- 
resentation which, again, relies on a mix of composite Gauss-Legendre 
quadrature and high-order analytic product rules 

While some steps in this scheme are standard, step 4, 6, and 7 are unique to 
the recently developed RCIP method. Conceptually, step 4 and 5 correspond 
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to applying a fast direct solver [26] locally to regions with troublesome ge- 
ometry and then applying a global iterative method. This gives us many of 
the advantages of fast direct methods, for example the ability to deal with 
certain classes of operators whose spectra make them unsuitable for iterative 
methods. In addition, this approach is typically much faster than using only 
a fast direct solver. 

Our numerical scheme can be extended to related problems of importance 
in e.g. band-gap structures, axially symmetric cavities for accelerators, and 
remote sensing of underground objects. Thus we can extend the method to 
scattering from homogeneous dielectric cylinders, scattering from multiple 
cylinders, scattering from cylinders in layered structures (c.f. [27J), scattering 
of plane waves at oblique angles from cylinders, and scattering from axially 
symmetric three-dimensional geometries. Some of these problems will be 
addressed in forthcoming papers. 
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